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model reduces to the conventional toric code in a uniform magnetic field. A 
quantitative analysis is performed for the perturbed Z3 toric code by applying 
a combination of high-order series expansions and variational techniques. We 
provide strong evidences for first- and second-order phase transitions between 
topologically-ordered and polarized phases. Most interestingly, our results also 
indicate the existence of topological multi-critical points in the phase diagram. 
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1. Introduction 

The concept of topologically-ordered quantum matter has been introduced by Wen in 
the context of high-temperature superconductivity [1, 2], and is crucial to characterize 
fractional quantum Hall states and topological insulators. The most striking property 
of topologically-ordered quantum phases is their dependence on non-local properties 
of the system. As a consequence, such phases cannot be characterized by a local order 
parameter, so that the celebrated Landau's symmetry-breaking theory cannot be used. 

More recently, topological order has attracted great interest in the field of 
quantum information due to its weak sensitivity to any local perturbation [3-.5]. 
Indeed, non-local degrees of freedom associated with this exotic order have been shown 
to be (topologically) protected against local sources of decoherence. This key idea is 
at the heart of topological quantum computation [3, (i]. It is thus of importance 
to quantify precisely this protection when perturbations are added. One prominent 
example where the effect of additional perturbations has been extensively discussed is 
the case of Kitaev's toric code [ ] . The toric code is an exactly solvable two-dimensional 
quantum spin model with a Z2 spin-liquid ground state possessing gapped (Abelian) 
anyonic excitations. It can be considered as one of the simplest models displaying 
topological order. Apart from the influence of temperature [7-11] and of disorder 
[12, 13], several works have investigated the effect of an external magnetic field in this 
model [14-20]. Interestingly, a very rich phase diagram containing first- and second- 
order phase transitions, multi-criticality, self-duality, and dimensional reduction has 
been found. Similar phase transitions out of topologically-ordered phases have been 
studied in the context of the Levin- Wen model [■M-23]. 

Actually, the toric code model can be defined for any discrete Abelian or non- 
Abelian group [3, 24, 25]. In this work, we present an extension of this model to 
Ztv degrees of freedom [20], which reduces to the conventional toric code for N — 2. 
Although excitations are still Abelian, qualitative and quantitative analyses of this 
model in the presence of local perturbations provide some insights in the understanding 
of topological phase transitions when more complex degrees of freedom are involved. 

The paper is organized as follows : in section 2, we describe an extension of 
Kitaev's toric code from Z2 to Zn degrees of freedom. Our starting point is Wen's 
plaquette model [27] for which a generalization to Z^r can be written down rather 
easily. The Z ^ plaquette model is then mapped onto a Z^r toric code which is shown 
to display topological order and anyonic statistics. 

To study the robustness of the topological order, we analyze the influence of local 
perturbations in section 3. For N — 2, such perturbations correspond to a uniform 
magnetic field. If the perturbations are strong enough, one expects conventional 
polarized phases that are not topologically ordered. As a consequence, a phase 
transition between the topological and the polarized phases must occur. Analyzing 
the breakdown of the topological phase is a very challenging problem for general N . 
Nevertheless, as we shall see for special kinds of perturbation, either exact mappings 
onto already known models exist or the model displays self-duality and dimensional 
reduction. 

In section 4, we focus on the case iV = 3 and probe the robustness of the 
Z3 topological phase for simple perturbations. To this end, we use perturbative 
continuous unitary transformations (pCUT) and a variational approach based 
on infinite projected entangled pair states (iPEPS). Note that this combined 
pCUT+iPEPS method has been already used successfully for the standard toric code 
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{N — 2) in an arbitrary magnetic field [j After a brief discussion of both methods 
and their combination, results are presented for two types of perturbation that can be 
viewed as natural generalizations of the N = 2 toric code either in a parallel field or 
in a transverse field. In the simplest case where the perturbation commutes with local 
charge (or flux) operators, we establish an exact mapping, valid at low-energies, onto a 
three-state clock model in a transverse field. For a ferromagnetic coupling, this model 
is known to display a weakly first-order transition (see for instance [2n]). However, the 
perturbation considered here also leads us to study the antiferromagnetic three-state 
clock model in a transverse field which, to our knowledge, has never been discussed 
in the literature. In this work, we find evidence for a second-order transition in this 
system (whose universality class still remains to be accurately determined). We then 
discuss the counterpart of an arbitrary parallel magnetic field in the N — 2 model. For 
such a perturbation, we obtain a rich phase diagram containing first- and second-order 
phase transition lines that form the boundary of the topological phase. Finally, we 
study the transverse-field problem that diplays dimensional reduction and self-duality, 
as the N = 2 model [IS]. Conclusions and perspectives are drawn in section 5. 



2. Construction of the generalized Zjv toric code 

In the following, we first introduce the Z^r-generalization of the plaquette model 
introduced by Wen [27]. The main reason to do so is that the plaquette model with 
Zjv degrees of freedom arises rather naturally from the Z2 conventional counterpart 
and it is very simple to write down. Afterwards we perform a mapping to a generalized 
toric code model with Z^r-Abelian anyons [26]. 



2.1. Zjv plaquette model 



We consider A^-state degrees of freedom located on the sites i of a square lattice whose 
unit-cell vectors ni and n2 are shown in figure 1 (left). The associated orthonormal 
states are denoted by \q)i, where g e Zat. Next, let us define the operators Zj and Xi 
as 



Zi\q)i = uj'^\q)i and Xi\q)i = \q-1) 



(1) 



where oj — . These unitary operators reduce to the conventional Pauli matrices erf 
and (7? for N = 2. On the same site, both operators obey the important "commutation 
relation" (Weyl algebra) 



XiZi — ujZiXi , 

which generalizes the well-known anticommutation relation cr|(7| — 
matrices. They obviously commute when acting on different sites. 
The Hamiltonian of the Z jv plaquette model is then defined by 



plaquette 



(2) 

-<jl(yf of Pauli 
(3) 



where Wp — 



ZbXb,ZtjXi^. Sites D, R, U, and L correspond to the four sites 
rii + n2,i + of an elementary plaquette p of the square lattice 



+ ni,i - 
(see figure 1). 

Using (1) and (2), it is easy to check that 



W W , 



p : p' 



= 0, 



(4) 
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Figure 1. Left : a piece of the square lattice on which one defines Wen's plaquette 
model. Right : a piece of the four-colored square lattice on which one defines 
Kitaev's toric code model with Zjv degrees of freedom. In Wen's model, degrees 
of freedom are located on the vertices of the lattice whereas in Kitaev's model, 
they live on the bonds of the lattice. 

for all p and p', so that all plaquette operators commute with the Hamiltonian. 
Furthermore, it is important to note that these operators obey = 1. In other 
words, Wp's are Z^v conserved quantities whose eigenvalues are simply {uj'',q G Zjv}- 
As for the standard Z2 plaquette model [27], this property ensures exact solvability 
of ifpiaquette- In the foUowing, we map the plaquette model onto a generalized toric 
code [3, 2G] as was already done for the Z2 case [29, SO]. The main advantage of the 
toric code is that the quantum statistics of the elementary excitations are simpler to 
identify. Additionally, it allows one to adopt an ad hoc language reminiscent of lattice 
gauge theories with Z^v degrees of freedom (see, e.g., [31, 32]). 




2.2. Zjv toric code 

2.2.1. Mapping In order to map model (3) onto a generalized toric code model [3], 
we introduce a translationally invariant and four-colored lattice as depicted on the 
right side of figure 1, where the A^-state degrees of freedom of the plaquette model are 
placed on the bonds of a square lattice. We then perform the following local, unitary 
transformations : 

V '> Z X > X 

iSbrown, magenta i ' iGycllow, green i ' 

brown, magenta i ' iGyellow, green i ' \ ) 

As a consequence, the plaquette operators Wp become different on the red and blue 
stars s and on the cyan and pink plaquettes p as illustrated in figure 1. To strengthen 
the analogy with the conventional toric code, we relabel the various operators as 
follows : 

%eblue ^ ^.eblue = XlX^XlX^ , W^;,p,^, ^ S,,pi„k = 44,44 • (6) 
Let us underline that the four-coloring of the lattice is mandatory if one wants to 
define Bp with Z (or Z'^) operators only. Nevertheless, other choices with smaller 
units cells are possible. 

Finally, we obtain the Hamiltonian of the Z^r toric code 
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Let us note that this model was aheady introduced in [ ] but differs from the Zjv 
toric code discussed in [3] which involves projectors Vs and Vp instead of Ag + Al and 
Bp + B^ (see next section for definitions). However, since both models are equivalent 
for N = 2,3, we shall (abusively) call them toric codes. 



2.2.2. Ground states and topological degeneracy A direct consequence of (4) is that 
all As and Bp operators commute with -ffxc- Thus, the ground-state energy (per site) 
Co is simply obtained by choosing the (possibly degenerate) minimal eigenvalue of the 
local operators —J{As + Al) or —J{Bp + B^), namely, Cq = —2Jcos{2nk/N). For 
J > and for any TV, the ground state is unique and obtained for k = 0. However, for 
J < 0, the ground state is unique for N even (in this case, one chooses k = N/2) but 
it is infinitely-many degenerate for N odd since, locally, one can choose k = (iV±l)/2. 
In the following, we will only consider the simplest case J > 0. 

Nevertheless, there are subtleties since the ground-state degeneracy also depends 
on the surface's topology as we shall now see on two simple examples. Let us first 
consider an infinite open plane for which no constraint on As's and i?p's exists. In 
this case, the ground state is unique and can be built as 



\gs)^Afl[Vsl[Vp\ref) , 



s p 

where A/" is a normalization constant and 

N 



1 



-1 



k=0 



N-1 



k=0 



(8) 



(9) 



Operators Vs (Vp) project on subspaces with eigenvalue 1 of the corresponding As 
{Bp). The reference state |ref) can be chosen arbitrarily provided it leads to |gs) 7^ 0. 

For instance, one may choose the fully-polarized state |ref) (x) |0)^ that already 



fulfills Vp |ref) — |ref) for all plaquettes p. 

Next, let us consider the Z^r toric code on a torus. In this case, there are two 
constraints 



1. 



II Bp 



1 . 



(10) 



so that the number of independent eigenvalues of As and Bp operators is reduced by 
two. However, as in the Z2 toric code, there exist conserved loop operators that can 
be chosen as : 



Z, = 



X, = 



n 



k iGCi , green 



kiGCs, green 




n 



iGCi .yellow 



n n 4 



Z2 = 



X2 — 



n 



k iGC2 , magenta 



n 



iGCs , yellow 



^i€C4 , magenta 




n 



z] 



iGC2 .brown 



n 4 



i$iC4 , brown 



where Ck with k E {1,2,3,4} are the non-contractible loops of the torus depicted in 
figure 2 (left). All these operators commute with i?TC, but only two of them can 
be chosen independently since [Zp,X^] ^ with /i S {1,2}. These two additional 
conserved quantities maintain the exact solvability of -ffxc • Concretely, if we choose 
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to label states with the eigenvalues — {w'} {q E Zn) of the operators Z^, we find 
that there exist N'^ ground states |gs, 21,^2) that can be written as 

|gs,..'Sa.'^) ^M'llVs {Xiy^ (4)''(8) |0>, , (12) 

s i 

where Af' is a normalization constant. More generally, following [3, 2()], one can show 
that for a compact surface with genus g, each eigenstate is TV^^-degenerate (at least) 
so that the system is indeed topologically ordered. 




Figure 2. Left : possible choices of non-contractible loops Ch with fc S {1, 2, 3, 4} 
for a system with periodic boundary conditions. Right : illustration of the semi- 
infinite strings Ss and Sp used to define charge and flux states in (13) and (14). An 
example of a counter-clockwise braiding contour for moving a charge qs initially 
located on a red star s around a flux qp located on a cyan plaquette p is represented 
in white (see text for explanations). The black circle locates the site (numbered 
8) where the crossing between the braiding contour and the string Sp of the flux 
occurs. 



2.2.3. Excitations and statistics Excitations of the toric code correspond to states 
which violate the condition that all eigenvalues of the As or the Bp operators are 
equal to 1. In other words, an elementary particle is a charge q on star s (a flux 
q on plaquette p) corresponding to an eigenvalue w"^, with q G and g 7^ 0, of 
As {Bp) (remember that q — defines the ground state). As As and Bp operators 
commute with i?TCi charges and fluxes are static excitations. Moreover, the form of 
the Hamiltonian implies that the energy of a many-particle state is simply the sum 
of the single-particle energies. In other words, charges and fluxes do not interact. 
In what follows, we give the explicit construction of single-particle states for general 
following the detailed construction given in [.30] for N — 2. As was already the 
case in the Z2 toric code, there is no local operator creating a single excitation. This 
is illustrated in figure 3 where one can see that an Aj (Zi) operator creates two 
excitations on neighboring plaquettes (stars). 

However, for an infinite plane with open boundary conditions, it is possible to 
consider single-particle states. Indeed, in such a system, a single excitation can be 
obtained by first creating a pair of charges (fluxes) and by taking one of the particles 
to infinity, at least in principle. It is clear that this is rather a gedankenexperiment that 
cannot be implemented practically (for instance, in a computer). Nevertheless, one 
may always consider a state where the two particles originating from the elementary 
pair-creation process are so distant that they are eventually independent. Thus, a 
one-particle state \q)^ with charge (flux) q on a ^ s (p) can be defined, for instance. 
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X 




as 



Figure 3. Illustration of the action of operators Xi (left) and Zi (right) on an 
eigenstate of As and Bp operators. A diamond on a star or a plaquette, with a 
it sign denotes a multiplicative change by ui^^ of the corresponding eigenvalue. 
The behavior depends on the site's color. 



\i)s= u zA\ n 4 igs), (13) 




\l)p= 11 11 lg^) (14) 

.magenta / brown / 

where the semi-infinite strings Sa are displayed in figure 2 (right). 

We shall now show on a specific but sufficiently general example that charges and 
fluxes obey mutual anyonic statistics. Let us consider an eigenstate, denoted by \ip), 
with a charge Qs and a flux qp at positions shown in figure 2 (right) . One can braid the 
charge around the flux along the counter-clockwise oriented white path drawn in this 
flgure, by acting on with the operator O = Zf^ . . . Z'^' (Z^"^ Zi" . . . Zf'. This can 
be checked from the action of Z^ operators shown in flgure 3 (right). Furthermore, 
from the definition of Bp operators given in (6), this operator O is nothing but the 
product of all operators for all plaquettes encircled by the braiding contour. 

Given that is an eigenstate of all plaquette operators with eigenvalue 1, except 
for the plaquette p where the flux is located and for which Bp \^) = ui'''' \tp), one gets 
O \ip) = 1-0). This non-trivial braiding phase is the signature of the mutual 

(Zn) anyonic statistics between charges and fluxes. It is similar to an Aharonov-Bohm 
phase, which explains the terminology of charges and fluxes employed to describe the 
excitations. The same argument allows one to show that a braiding of a charge (flux) 
around a charge (flux) leads to a trivial phase which is reminiscent from the bosonic 
statistics of charges (fluxes). In addition, it is clear that there is a hard-core constraint 
since one cannot create two particles on the same star or on the same plaquette. 

Let us end this discussion by showing that the phase can be obtained in another, 
complementary way. One can write explicitely the state {ip) as 




1^)- n n ^1 n n 

.green / WtEtSs .yellow / .magenta / brown / 

Then, one can compute O by commuting O with the operators appearing in the 
above expression using (2) and O |gs) = |gs) since the ground-state is flux-free. For the 
particular braiding represented in flgure 2 (right), the non-trivial phase will appear 
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from the commutation Z^^'Xl'' — ui '^"'^p X'^" Zf" at site number 8 where the braiding 
contour and the string of the flux Sp intersect. 

3. Ztv toric code in the presence of local uniform perturbations 

In this section, we first define the general local perturbations for the Zat toric code 
model. The presence of such perturbations destroys the exact solvability of the toric 
code model. As mentioned in the introduction, a phase transition has to take place 
when the perturbation strength increases since, for J = 0, the ground state is fully 
polarized and thus not topologically ordered. Thereafter, we focus on special examples 
that allow for a detailed investigation of this transition for general TV. 

3.1. General structure of the perturbation 

Here, we shall only consider perturbations which act locally on a site i. A basis of 
the space of local unitary operators can be conveniently written down in terms of the 
operators Xi and Zi and their powers. Hermitian combinations of these operators 
lead to the following general form of local uniform perturbations 

H„n = -Y. (^M»^i^r + h-c.) , (16) 

i 

where hi^„i G C, and {l,m) e Z^ with ^ ^ ^ N/2. The action of the operator 
X\Z™ is displayed in figure 4. As can be inferred from this figure, the perturbation 
violates the local conservation of charges (fluxes) when m ^ ^ 0) since it induces 
creation and annihilation of pairs of excitations as well as hopping processes. However, 
it violates neither the conservation of the total charge Ylis 9s modulo TV, nor the 
conservation of the total flux qp modulo N [we recall that a star s (a plaquette p) 
is said to carry a charge qs (a flux qp) if As {Bp) has eigenvalue w'" (w''''), with qs and 
qp being defined modulo N ] 




Figure 4. Illustration of the action of the operator X\ZV^ on an eigenstate of Ag 
and Bp operators. A diamond on a star or a plaquette, with a value k = itZ, itm, 
denotes a multiplicative change by oj* of the corresponding eigenvalue. The 
behavior depends on the site's color. 

In fact, these conservation rules can be more conveniently rewritten as 
conservation rules of "unphysical" total charge and flux, belonging to Z and not 
to Ztv (they will also prove to be useful later on, see section 3.3). Let us 
denote the "unphysical" charge (flux) at star s (plaquette p) as qs {%)■ They are 
related to the true/physical charge (flux) via the following relations qs = qs mod N 
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ilp — Qp mod N). They can be defined non ambiguously by "fixing a gauge" as 
we explain now (having in mind a purely computational perspective). Given one 
eigenstate of Htc, it is possible to choose any value for the "unphysical" charges 
and fluxes (provided they yield the correct physical charges and fluxes). Then, the 
full Hamiltonian Htc + X); m Hi,m can be studied in the subspace of eigenstates of 
Htc, spanned by the repeated action of the perturbation. Each of these states can be 
assigned unique values of qs and qp using rules shown in figure 4. It should be clear 
that Xg qs and qp are the same for all states in the generated subspace, which 
means the Hamiltonian is block diagonal and the total charge X^s 5s and total flux 
Xp 9p are conserved integers. 

3.2. Simplest cases 

Let us first discuss the case I = for which the perturbation acts only non trivially on 
charges (which arc no longer locally conserved) whereas fiuxes remain static gapped 
excitations that might be seen as sources of Aharonov-Bohm-likc phases for the moving 
charges. Of course, the respective roles of fluxes and charges are exchanged if / ^ 
and m = 0. 

In this work, we are interested in transitions between the topological phase 
(existing for small enough perturbations) and the polarized phase expected for 
J = 0. To get a flrst idea about the associated physics, let us consider the simplest 
perturbation corresponding to {I = 0,m = 1), for which the Hamiltonian of the 
perturbed toric code reads 

Htc + Ho,i = - (^^ + - (^p + ^ ^^Y. {^i + ^l) ' (1^) 

s p i 

where we set /iq.i = /iz G M . Since we are interested in the low-energy properties 
and since [Ho,i,Bp] = for all p's, we only consider the flux-free subspace where all 
-Bp's have eigenvalue 1, in which the energy of the fluxes is minimal. Then, following 
the procedure discussed in [14, 15] for the special case N = 2, let us denote by \q)s 
the eigenstates of Ag's with eigenvalues w'. From figure 3 (right), one can check that 
the action of + -^1 on a site i located between two neighboring stars s and s' 
is equivalent to X^X^, + XjX^, where, following (1), we introduce the operator Xg 
defined by Xs\q)s = \q~ l)s- Thus, defining Zg — Ag (such that Zs\q)s = uj'^\q)s), one 
can map Hamiltonian (17) onto the A^-state clock model in a transverse field [-i'd] 

i?ciock = -2JNp^jY,{Zs + Zl)-hz (X.XI + XtX,,) , (18) 

s {s,s') 

where (s,s') denotes nearest-neighbor stars s and s'. The term —2JNp arises from 
the replacement of all Bp-operators by their eigenvalue 1 (Np denotes the total 
number of plaquettes). It is important to stress that this mapping preserves neither 
the degeneracies of the energy levels (hence the topological order) nor the quantum 
statistics. However, the zero-temperature phase diagrams of Htc + Hq i and iJciock 
are exactly the same. 

Let us remark that the coupling term hz in (18) stems from the local perturbation 
in (17) that can be either positive or negative. When hz > 0, the coupling between 
stars in -ffciock is ferromagnetic whereas hz < leads to antiferromagnetic interactions. 
This distinction is irrelevant for even N since, in this case (and for a bipartite 
lattice), one can always perform local unitary transformations which map -ffciock onto 
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a ferromagnetic model. By contrast, for odd N, one must distinguish between both 
signs that may lead to various types of transitions. 

Unfortunately, few results are available in the literature concerning the two- 
dimensional quantum clock model in a transverse field except for N = 2 (Ising model) 
where a second-order transition occurs, for = 3 (Potts model) where a weakly first- 
order is expected for hz > (see for instance [ ]), and for = 4 where the model 
is equivalent to two decoupled transverse-field Ising model [23]. In such a context, 
the second-order transition found in section 4 for = 3 and hz < opens some 
interesting perspectives. 

3.3. Self-duality 

The Z2 toric code in a transverse field is known to be self-dual [IS, 34]. We shall 
now show that this property still holds for a general value of N, provided one chooses 
particular values of I and m. Let us consider the Hamiltonian H^c + -ff;.m [see (7) 
and (16)]. This Hamiltonian will be self-dual if its spectrum is symmetric (up to 
degeneracies) under the exchange J O [/i/^mj- Roughly speaking, this will be ensured 
provided the "roles" of Ag and Bp operators can be played by the operators X^Z™, 
once stars and plaquettes are exchanged with sites, i.e., when considering the dual 
lattice illustrated in figure 5. 

In the limiting cases, J — and hi^m — 0, self-duality imposes that spectra of i?TC 
and Hi^rn are the same (up to degeneracies). Setting /i;.,„ — Ift./.mje"^' '" , this means 
that e"^' "AT^Z™ must have the same spectrum as Ag (or Bp). Since this spectrum 

is {uji,q e Zn}, this leads to two constraints : (i) {e"f''-^X\ZV^)^ = 1 and (u) 

(e'-^'-'X^Zr")" 7^ 1, V n e {1, . . . , A^-1}. Noting that (X^Zf)" = a;''"'^^^^ XfZf", 
the first constraint reads 

e'^-^'^-w'"^-^ = 1 , (19) 

whereas the second one can be rephrased as 

7^ mod iV or mn 7^ mod TV, Vn G {1, . . . , - 1} . (20) 

The third and most stringent condition is that the action of Hi m on eigenstates 
of TJtCj of all Ag and Bp operators (see figure 4), is the same as the action of 
i?TC on eigenstates of ff;,™, i-c, of all X\ZV^ operators (see figure 6). After noticing 
that the numbers of minus signs, as well as their exact positions, can be made the 
same in both figures by exchanging the roles of X\Z™ and its hermitian conjugate on 
magenta and green sites (or on yellow and brown sites), this last condition reads 

m = l or m = N — I. (21) 

For N = 2, one recovers the model studied in [1^] (up to a factor of 2). Indeed, 
the only solution to the above equations is Z = m = 1 and = ±7r/2. Thus, the 
perturbation reads ±2hy cr?', while the toric code Hamiltonian can be simplified to 
—2J{J2s ^s + J2p ^p) since the star and plaquctte operators are Hermitian for N = 2. 
Note that the present analysis also shows that the sign of the field term is irrelevant. 

For the case A^ = 3, to which the whole next section is devoted, condition (19) is 
fulfilled for e^"^' " = l,w and a;^, whereas (20) and (21) are fulfilled for (^ = l,m = 1) 
and (? = l,m = 2). The six corresponding perturbations yield the same spectrum so 
that, for simplicity, we shall only consider Hi i with (j>i i = 0. 

To conclude this discussion, let us show that self-duality in these models is also 
responsible for additional symmetries, as was already observed for the special case 
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N = 2 [18]. Consider figure 4 with m — I (the same discussion holds for m — N — I). 
It is clear that the sum of the numbers appearing in the diamonds, on diagonals or 
anti-diagonals, is either 0, 21 or —21, which is always even. This allows one to define 
conserved parity operators such as for example (— l)2Isercd fs+I]p(=pi„k ip ^ where the sum 
runs overs red stars and pink plaquettes forming a given diagonal. In order for this 
parity operator to be conserved when is odd, one has to consider the charge and 
flux numbers Qs and qp belonging to Z, introduced at the end of section 3.1 instead 
of Qs and Qp that belong to Z^r. Of course, a dual discussion can be given by working 
on the dual lattice and in the eigenbasis of X\Z™ operators. As in the Z2 toric 
code in a transverse field, the conservation of these parity operators constrains the 
dynamics, and the model displays dimensional reduction. This dimensional reduction 
was originally discussed in the Xu-Moore model [-35-37] which has the same spectrum 
as the Z2 toric code in a transverse field. A detailed discussion about these issues can 
be found in [38]. 




Figure 5. Transformation of the original lattice (left) to the dual lattice (right). 
Plaquettes and stars become sites, and vice versa. 




^sGblue 



Figure 6. Illustration (on the dual lattice defined in the right part of figure 5 
) of the action of the operators As and Bp on an eigenstate of X^ZV^ operators 
for all sites i. As in figure 4, a diamond on a star or a plaquette, with a value 
k = ±1, ±m denotes a multiplicative change by oj* of the corresponding eigenvalue. 
The behavior depends on the site's color (so on the star and plaquette types on 
the original lattice). 
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4. Perturbing the Z3 toric code 

4.1. Model 

In this section, we focus on the case = 3 and we study the robustness of the Z3 toric 
code with respect to simple perturbations. To this aim, and for the sake of simplicity, 
we consider the following Hamiltonian 

H{hx,h^,hz) = Htc + Hi,o + Hi^i + Ho,i, (22) 

s p i 

-h^Y. {^i^i + ^1^1) -^^T. {^i + ^1) ' (23) 

i i 

where we choose J = 1/3 in order to set the elementary excitation gap of the 
unperturbed Hamiltonian H^c to unity. In addition, we restrict our discussion to real 
parameters and set /ii.o = hx, = h±, ft,o,i = hz- For N = 2, this Hamiltonian 
(with the proper phase factor for Hi i discussed in section 3.3) corresponds to the Z2 
toric code in a uniform magnetic field studied in [ ; ] so that this choice is well suited 
for a comparison between both systems. In the following, we adopt a language similar 
to that used for N = 2. Thus, Hi^q and Hq i (-ffi.i) will be considered as "parallel" 
("transverse") perturbations. 

The determination of the full three-dimensional phase diagram of Hamiltonian 
H{hx,h±,hz) is a difficult problem. Therefore, as was done initially for the Z2 case 
[17, 18], we will study parallel and transverse cases separately, but let us first discuss 
the methods used. 

4.2. Methods 

As already discussed, the system has to undergo phase transitions when perturbed 
with local operators. As in the perturbed Z2 toric code, one expects first- and second- 
order transitions in the phase diagram [ ] . In order to analyze the breakdown of the 
topological phase, we combine the pCUT method and the variational iPEPS algorithm. 
This approach is motivated by the fact that the pCUT gives reliable estimates for 
second-order phase transitions whereas the iPEPS algorithm, as a variational tool, is 
especially sensitive to first-order transitions. 

4-2.1. pCUT The method of continuous unitary transformations has been introduced 
in reference [39] and general aspects of its perturbative variant pCUT can be found in 
reference [40]. In what follows, we focus on points that are specific to the application 
of the pCUT method to a topologically-ordered phase [17-1!)]. 

To apply the pCUT, it is essential that the unperturbed Hamiltonian (here i?Tc) 
possesses an equidistant spectrum ^ that is bounded from below [40]. These two 
constraints are satisfied in the Z3 toric code as long as gaps of charges and fluxes are 
identical. In this case, one can interpret the toric code as a counting operator Q of 
charges and fluxes in the system 

Htc= - I {Ns + Np) + Q, (24) 

X We would like to stress that the spectrum of the Zjv model studied in this paper is only equidistant 
for N = 2,3,4, so that one cannot use this machinery for other values of A'^. However, it could be 
applied to Kitaev's Zjv toric code [■>], whose spectrum is equidistant for all N. 
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where Ng {Np) denotes the total number of stars (plaquettes). Therefore the constant 
term represents the ground-state energy (remember that we set J = 1/3). 

It is then possible to rewrite the local perturbations as J2n where r„ changes 
the particle number in the system by n, i.e., [Q,T„] = nTn- The pCUT maps, order 
by order in the perturbation, the Hamiltonian H onto an effective Hamiltonian 
unitarily equivalent to H (same spectrum but different eigenstates) , that reads as 
follows in the eigenbasis of the bare Hamiltonian i/TC • 

2 °° 
H-^ = --{Ns + Np) + Q + C(mi,...,mfc)r,„, ...T™, . (25) 

k — l mi-f...+mfc— 

As explained in [ i' )], the coefficients C(mi, . . . , rn^) are model-independent rational 
numbers. An essential property of the effective Hamiltonian is that \H'^^ , Q\ = 0. As 
a consequence, the number of quasiparticles (QPs) in the system, i.e., eigenstates of 
Q, is a good quantum number. In the perturbed toric code, QPs are dressed anyons 
adiabatically connected to the corresponding bare charges and fluxes. 

To determine the zero-temperature phase diagram, we focus on the low-energy 
spectrum of . Essentially, one must study in the 0-QP subspace (to 

compute the ground-state energy) and in the 1-QP subspace (to obtain the low-energy 
gap A). We emphasize that computing the low-energy gap from the 1-QP subspace 
is meaningful as long as there are no bound states with lower energies. This is a 
working hypothesis that is crucial in what follows. As discussed in section 2.2.2, for 
open boundary conditions, the ground-state is non degenerate and the ground-state 
energy is nothing but the expectation value of H^^ on the ground state of the bare 
Hamiltonian 

= (gs| ff'^ff |gs) . (26) 

The structure of the 1-QP subspace is less trivial. Indeed, for iV = 3 there are four 
different kinds of excitations : charges (g^ = 1) and anti-charges {qs — —1 — 2 mod 3) 
living on stars, as well as fluxes [qp — 1) and anti-fluxes [qp = —1 = 2 mod 3) living 
on plaquettes. The associated subspaces are not connected by because of the 
symmetry ensuring that the total charge and total flux are conserved (modulo = 3 
when working with the physical charge and flux). Thus is a 1-QP hopping 

Hamiltonian in each of these sectors, and is therefore easily diagonalized, once the 
hopping amplitudes are determined. One then obtains dispersion relations for the 
four kinds of excitations, which only give two different energies because charges and 
anti-charges, as well as fluxes and anti-fluxes, play symmetric roles. Finally, one can 
compute the gap A as the minimum over all momenta of these two energy bands. 

Let us point out that the major challenge, in both the 0-QP and 1-QP sectors, lies 
in the computation of matrix elements of Indeed, one has to take into account 

the non-trivial braiding phases coming from virtual fluctuations of the excitations. 
Although the method yields results in the thermodynamical limit, the linked-cluster 
theorem (see for example [41, 42] and references therein) allows one to compute the 
hopping amplitudes by considering flnite-size clusters (albeit one needs a growing 
number of them when the order of the perturbation increases). 

At the end of the day, one obtains a high-order series expansion of the ground- 
state energy Eq and of the 1-QP gap A. The extrapolation of A with standard 
resummation techniques (see e.g. [43]) allows a reliable determination of second-order 
phase transition points and thus of the boundaries of the topological phase (assuming 
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that bound states of elementary QPs are not relevant and do not have an energy 
smaller than that of a single QP). 

Unfortunately, as already stated, series expansions are not adapted to detect 
possible first-order phase transitions (in particular when having series in one phase 
only) so that one needs a complementary tool which we now describe. 



4-2.2. iPEPS The so-called iPEPS algorithm [ ] is a variational method which, 
as such, is aimed at approximating the ground state of two-dimensional quantum 
lattice systems by employing a tensor-network approach. Details about this method 
have already been extensively discussed in the literature [44-48]. For completeness, 
though, we explain some of the basic features of the algorithm, focusing on those that 
are relevant for the study of the perturbed Z3 toric code. 

In the iPEPS algorithm, the quantum state of the infinite square lattice 
is represented by a projected entangled pair state (PEPS) [44, 49]. In the present 
problem, we choose a PEPS with four tensors denoted P, Q, R, and S per unit cell 
(see figure 7). Each of these tensors depends on 0{dD'^) complex coefficients, where 
c? = iV = 3 is the dimension of the local Hilbert space at each site, and D is the so-called 
bond dimension of the PEPS. This bond dimension controls the maximum amount of 
entanglement carried by the PEPS wave function and, consequently, the accuracy of 
the ansatz. Following the discussion in [ ] for N = 2, one can show that the ground 
state of the (non-perturbed) Zn toric code is a D = A^-PEPS. Obviously, for J = 
the fully-polarized ground state is also a (trivial) PEPS {D — 1). Consequently, at 
least in these two limiting cases, PEPS are exact ground states of the Hamiltonian. 

In practice, for a given D, the goal is to find the coefficients of tensors P, Q, R and 
S that best approximate the ground state of H. These coefhcients can be determined 
by an imaginary-time evolution driven by the Hamiltonian, since 



where l^'gs) is the ground state of H and |^'o) is any initial state that has a non- 
vanishing overlap with the ground state. The approximation of this evolution is 
performed in a similar way as explained, for instance, in [44, 47]. 

(i) The whole evolution is splitted into small imaginary-time steps 5t by using a 
Suzuki- Trotter expansion of the evolution operator e~'^^ . More precisely, writing 
the Hamiltonian as a sum of four-body terms 

H = Y, /i''^'"'' > (28) 

ijkl 

we consider, at each time step, the action of the fonr-site operator 

^[i,-fcJ]^g-6rhl--l^ (29) 

(ii) At each imaginary-time step the state is approximated by some PEPS with 
the considered structure and bond dimension D. For instance, if at step r we 
have a PEPS |*(r)), then the evolved state |«'(t + St)) = |5'(r)) is also 
approximated by a new PEPS |v1/(t -|- St)) with the same structure. Practically, 
this approximation is achieved by minimizing the distance |||5'(T-|-(5r)) — |^'(T-f 
5r))||^ with respect to the coefficients of tensors P, Q, R and S of the new PEPS. 
In our case, we have carried out this minimization simultaneously over the four 
tensors by using a standard conjugate-gradient algorithm. 




Figure 7. Schematic representation of the (translation-invariant) PEPS 
considered in this study with four different tensors [P, Q, R, and 5) per unit 
cell. Tensors are represented by circles, and their indices by lines. Lines that 
connect different circles correspond to bond indices shared by two tensors and 
can take up to D different values. Open lines correspond to physical indices, 
which take d = N = 3 values (dimension of the local Hilbert space on each site). 

Quite importantly, step (ii) as well as the evaluation of expectation values of local 
observables involves the contraction of an infinite two-dimensional tensor network. 
This contraction can be approximated by various schemes [44, 47, 51]. Here, we 
choose the Directional Corner Transfer Matrix Approach introduced in [47] that can 
be easily adapted to deal with different types of two-dimensional tensor networks [4(S] , 
including those considered here. An important parameter in these manipulations is 
the so-called bond dimension of the environment x which controls the accuracy of the 
approximations involved at this step. 

Thus, according to the discussion above, there are four possible sources of error 
in the iPEPS algorithm. 

(i) The size and the shape of the considered unit cell. The error is reduced when the 
unit cell gets larger. 

(ii) The bond dimension D of the PEPS. The error is reduced when D gets larger. 

(iii) The imaginary-time step St. The error is reduced when 5t gets smaller. 

(iv) The bond dimension of the environment x- The error is reduced when x g^ts 
larger. 

In this study, we have fixed D = 3, x — 30, St = 10^^, and the aforementioned 
"four-site" unit cell. We checked that the increase of precision obtained by varying 
the values of the parameters is within the error bars obtained by the pCUT approach. 
Therefore, this choice of parameters turns out to be sufficient for our purposes. 

4-2.3. pCUT+iPEPS To determine the boundaries of the topological phase, we 
combine results from pCUT and iPEPS algorithm as explained below. To simplify 
the discussion, let us assume that we have only one control parameter h> Q. 

Let us recall that a second-order transition is associated to the closure of the 
low-energy gap. Assuming that the relevant gap comes from the 1-QP sector, the 
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critical point h"^ is then defined by A {h'^) = 0. As already mentioned, this point 
can be efficiently computed with the pCUT method by extrapolating the high-order 
series expansion of A. Typically, with the maximum orders reached in this study, one 
can determine /i'^ with a relative precision of the order 10~^-10~'^. However, if some 
level crossing occurs, one faces a first-order transition that cannot be captured by the 
criterion A = 0. This is why it is crucial to use the iPEPS algorithm to compute, 
variationally, the ground-state energy. Indeed, denoting by Bq^^'^ and ejf^^^ the 
ground-state energies calculated by both methods and assuming the existence of a 
point h* where ejf^^^ < eg'^'^"'", two cases must be considered. If h* > h'^, it means 
that the iPEPS algorithm does not detect any level crossing for the ground state before 
the critical point, and hence a second order transition is likely taking place at /i'^. By 
contrast, if h* < /i'^, then it means that a level crossing definitely occurs before the 
critical point h'^ and we conclude that a first-order transition takes place at h* . 

Of course, this reasoning would be exact if (i) one would have an infinite precision 
in both methods which is clearly not the case and (ii) no bound-state with lower 
energy exists (see discussion above). The accuracy of this combined pCUT-t-iPEPS 
approach is limited by several sources. First, the series expansion is performed up 
to a finite maximum order and the error of resummation schemes like the dlog-Pade 
extrapolation is hard to quantify. Second, the variational iPEPS algorithm is limited 
by the values of the parameters D, x, St as well as the structure of the tensor network 
itself (see section 4.2.2). Additionally, it is numerically challenging to extract the 
global minimum of the variational ground-state energy. Reasonably, one can state 
that the combined pCUT-|-iPEPS approach works well as long as the error bars of 
both methods when determining the values h* and h'^ are small compared to the 
difference \h* — It!^]. As explained above, in the present work, this relative error on h* 
and h'^ is of the order 10^^-10"^. 

4-3. Results 

Let us now present our results concerning the perturbed Z3 toric code H{hx; h±, hz) 
[see (22)] for several limiting cases. 

4-3.1. The case hx — h± — As already mentioned in section 3.2, the low-energy 
physics of the perturbed toric code H{0, 0, hz) corresponds to the iV-state clock model 
in a transverse field (18). For = 3, this model is also equivalent to the three-state 
Potts model in a transverse field. The extension of the topological phase can therefore 
be obtained by directly analyzing this model, which is simpler since one only has to 
consider charge degrees of freedom that live on stars of the square lattice. Indeed, 
remind that fluxes are frozen (conserved) and even absent in the low-energy sector. 

As a first step, apart from the pCUT-|-iPEPS analysis, let us perform a standard 
mean-field calculation that, as we shall see, already captures the main qualitative 
aspects of the phase diagram although it relies on a trivial (non-entangled) variational 
state. To this aim, let us consider the following trial wave function 

l*)= <S> l^b), (g) l^r),, with l^b,), -ab,.|0),+6b,r|l),+Cb,r|2), , (30) 
s^Ebluc sGrcd 

where the coefficients ab,r, 6b, r, and Cb,r are chosen such that the local wave 
functions \tph,r) ^ are normalized and minimize The introduction of different 

wave functions on red (r) and blue (b) stars is needed to accomodate with both 
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ferromagnetic order and anti-ferromagnetic order (expected for hz > and hz < 
respectively). For the clock model, the sublattice magnetization 

Mb,r (hz) - ((^ba-l Xs IVb.r) (^b,r| Xj |Vb,r)) ' , (31) 

is a proper order parameter. The topological phase of the perturbed Z3 toric code 
corresponds to the disordered (symmetric) phase of the three-state clock model 
characterized by Mb,r = 0. This phase is obtained for \hz\ ^ J = 1/3. By 
contrast, an ordered (broken) phase with Afb,r 7^ is expected for large perturbations 
\hz\ ^ J — 1/3. Once again, let us emphasize that the sign of hz is important for 
N = 3. Indeed, for large positive hz, a polarized phase (with all stars in the same 
state) is stabilized, whereas for large negative hz a "staggered" order arises. 

Our results for the mean-field order parameter are summarized in the upper panel 
of figure 8. For hz > 0, we find a first-order transition at /i^^^ = 1/9. At this point, 
the order parameter jumps discontinuously from to 1/2. For hz > ft.^'^^, the system 
is uniformly polarized (M^ = Afb) and lim Mb^r = li a-s it should. For hz < 0, we 

obtain a continuous quantum phase transition at l?^^ = —1/8. As discussed above, 
the ordered phase found for hz < h^ < spontaneously breaks the translation 
symmetry of the system, i.e., A/b 7^ Afj. In the mean-field aggroximation, both 
(sublattice) order parameters and vanish as (/i^^^ — hz)^ with /3^^ = 1/2. 
Furthermore, as can be shown by a simple first-order perturbation theory in J/hz, 
one has lim Afb = 1 and lim = 1/2 (or vice versa). 

To go beyond this mean-field analysis and to obtain more quantitative results, let 
us now turn to the pCUT-|-iPEPS analysis. As the Hamiltonian (18) only contains 
one-site and two-site terms, the evaluation of effective matrix elements in the pCUT 
calculation can be performed using a full graph decomposition [43, 52] that allowed 
us to reach order 11 (see Appendix A). The iPEPS algorithm also considerably 
benefits from the absence of four-site terms in i/ciock- As can be seen in Fig. 8, 
the PCUT-|-iPEPS method is in qualitative agreement with the mean-field treatment 
since it also predicts a first-order transition for hz > and a second-order transition 
for hz < 0. However, for hz > where our results match those given in !] for the 
Potts model (up to a rescaling), we find h*^ ~ 0.126 and /i^ — 0.129. Thus, within the 
combined pCUT-|-iPEPS scheme, we are led to conclude that a (weakly) first-order 
phase transition occurs at h*^ ~ 0.126 in agreement with the results given in [28, 4()]. 
As can be seen in Fig. 8 (upper panel), the calculation of the order parameter using the 
iPEPS algorithm confirms a first-order behavior (jump of the magnetization). Note 
also that the relative difference between /i^ and h^^^^ is about 12%. 

The situation is different for hz < 0. In this case, we find h* ~ —0.204 and 
ft,^ — —0.195 (see table Al for more details). Let us stress that we observed a very 
good agreement between iPEPS and PCUT calculations with a relative error between 
both ground-state energies smaller than 10"'' for h% < hz < 0. We therefore conclude 
that, within our scheme, a second-order phase transition occurs at /i^ — —0.195. In 
this parameter region, the discrepancy between /i^ and /i^'^^ = —1/8 is larger than 
35%. To the best of our knowledge, this second-order phase transition has never been 
discussed in the literature and it is therefore very desirable to further characterize 
its universality class. Unfortunately, the iPEPS calculation of the critical exponent /? 
associated to the order parameter is known to be specially sensitive to finite-D effects. 
As shown in the one-dimensional Ising model in a transverse field [54], this exponent 




Figure 8. Comparison of pCUT, iPEPS, and mean-field results for the N = 3 
clock (Potts) model in a transverse field hz {J = 1/3). At order 11, all dlog-Pade 
approximants agree within the width of the black lines. Vertical dashed lines 
indicate the position of defined by A(/i^) = 0. Vertical dotted lines indicate 
the position of /i^ beyond which Cq'^^^ < e'^'^^'^ . Upper panel : sublattice 
magnetizations Mb ^ as a function of hz- Central panel : 1-QP gap A as a 
function of hz- Lower panel : ground-state energy per site eq as a function of 
hz- Inset : zoom of eo close to the critical point for negative hz- 



also eventually reaches its mean- field value 1/2 for any finite D when approaching 
the critical point. Nevertheless, the exact value (1/8 in the latter problem) can be 
observed in a field range near the critical point whose size increases with D. In a 
two-dimensional system, it is very difficult to increase D in order to perform a similar 
study. For the problem at hand, our fixed bond dimension D = 3 only allowed us to 
see an exponent different from /3^^ = 1/2 in a very small region, and this was not 
sufficient to determine a reliable value. 

Alternatively, dlog-Pade extrapolations of high-order series expansion of A allows 
one to determine the exponent driving the closure of the gap at the critical point. More 
precisely, denoting z the dynamical exponent and v the correlation length exponent, 
one has A ~ |/i — /I'^j^" for h near h"^ (see e.g. ]). At order 11, we found z ~ 0.71 
but, as can be seen in table Al, it is clear that this value is poorly converged. However, 
to roughly estimate the error, one may draw a parallel with the N = 2 problem (i.e., 
the Ising model in a transverse field) for which series expansion of the gap have been 
computed up to order 13 in [ ]. There, using the order 11 results, one gets an 
exponent z ~ 0.645 which only differs by a few percent from the commonly accepted 
values z = 1 and v = 0.630(1). Hopefully, for A'' = 3, we are also close from the 
asymptotic value but one clearly needs a more quantitative study to clarify the nature 
of this quantum phase transition. From that respect, Monte-Carlo simulations could 
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provide valuable insights. 

4-3.2. The case h± = Let us now turn to the case when both perturbations i?o,i 
and Hi Q are present so that neither charges nor fluxes are locally conserved anymore. 
Thus, one has to treat charges and fluxes at the same level and to carefully take 
their mutual statistics into account in the virtual braiding processes. This constraint 
strongly reduces the maximum order that can be reached in the pCUT calculation. 
Instead of order 11, we only computed eo and A up to order 7 for this case (see 
Appendix A). The iPEPS calculation becomes more involved as well since one now 
has to deal with 4-body interactions (instead of 2-body interactions in the three-state 
clock model). In other words, our results are less accurate when a more complex 
perturbation is considered, as expected. 

In figure 9 (left), we display the phase diagram obtained by combining pCUT 
and iPEPS algorithm following the same procedure as previously but for arbitrary 
directions in the {hx,hz) plane. As can be seen, the shape of the topological phase 
is symmetric under the exchange of hx hz- This is due to the fact that the 
perturbation -ffo.i + -^1,0 respects the "charge-flux symmetry" present in H^c- The 
transition lines that mark the boundaries of this phase are found to be either first-order 
or second-order lines. 

Let us start with a discussion of the (two symmetric) first-order (cyan) lines that 
are directly connected to the weakly first-order transition points (cyan diamonds) of 
the three-state clock model in a (positive) transverse field. Near these points, h"^ 
and h* are found to be very close so that it is challenging to clearly decide whether 
the transition is first- or second-order. However, given the precision reached and the 
existence of some well-identified first-order points (not shown in figure 9) away from 
the cyan diamonds, it seems reasonable to argue that the two lines are likely first 
order. 

The situation is different for the other part of the topological phase whose 
boundaries are connected to the second-order points (black diamonds) of the three- 
state clock model in a (negative) transverse field. Here, the pCUT+iPEPS approach 
is always clearly consistent with a second-order phase transition. 

In this region, the intersection of the two second-order transition lines (green 
triangle) is reminiscent of the phase diagram of the perturbed Z2 toric code in a parallel 
magnetic field where a multi-critical point was found [IG, 17]. In the latter case, this 
crossing point is also connected to a finite-length first-order line that lies outside the 
topological phase. For the present problem {N = 3), we have not performed a small- J 
perturbation theory likely to reveal a similar feature. Nevertheless, we computed the 
gap exponent z v along the second-order (black) lines. As can be seen in figure 9 
(right), the behavior of this exponent is rather flat (same values as for hx = or 
hz = 0) except at its extremities where it increases signiflcantly. This may indicate a 
different universality class at the crossing points (green triangle and magenta squares) 
as also found in the N — 2 problem. Anyway, let us stress that there are some 
limitations of our approach for the current problem. First, the perturbative expansion 
at order 7 is not sufficient to determine the gap exponent z v accurately. Second, 
the finite width in h^ where the gap exponent differs from the value at h\ = is 
likely an artifact of the finite-order series. Indeed, we rather expect that all points 
except crossing points belong to the same universality class as the anti-ferromagnetic 
clock model in a transverse field but one cannot exclude other scenarios. Once again, 
to obtain more quantitative results, it would be very valuable to perform numerical 
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Figure 9. Left : zero-temperature phase diagram of H{hx,0,hz) obtained 
by the pCUT-|-iPEPS approach. Cyan (black) lines denote first-order (second- 
order) phase transitions. The first-order (second-order) phase transitions for the 
ferromagnetic (anti-ferromagnetic) three-state clock model in a transverse field are 
marked by cyan (black) diamonds. Magenta squares locate the crossings between 
first- and second-order transition lines while the crossing of the two second-order 
lines is taking place at the green triangle. Right : plot of the gap exponent zu a,s 
a function of /i^ along the horizontal second-order phase transition line (see left 
part). The black solid line is the result obtained from the order seven expression 
of the 1-QP gap. The black square is our best estimate of 2 ^ for the simpler case 
of one parallel perturbation using the order eleven series expansion obtained for 
the three-state Potts model in a transverse field (see section 1.3.1). 



simulations of tliis model by means of alternative methods. 

4-3.3. The case hx = hz = To conclude this analysis, let us consider the case of 
a "transverse" perturbation that, as discussed in section 3.3 for general N, leads to 
a self-dual spectrum for H{0, h±,0). Since J > and N is odd, one must also have 
h± > 0. A convenient parametrization for this problem consists in setting J = ^ cos 
and h± — ^ sin 9 with 9 E [0, 7r/2]. The unperturbed toric code corresponds to 9 — 
whereas for 9 = tt/2 the Hamiltonian is purely local. The self-dual point is located 
at 9 = 7r/4 (J = h±) so that the spectrum is symmetric under the transformation 
9 O tt/2 — 9. Thanks to self-duality, one can determine directly the nature of the 
transition by studying the singularity of the ground-state energy. Note that we could 
not use this criterion in previous cases since we did not have reliable informations 
outside the topological phase. As can be clearly seen in figure 10 (lower panel), the 
ground-state energy displays a kink at the self-dual point (^^|e=7r/4 is discontinuous) 
so that a first-order transition occurs there. 

As explained in section 3.3, the parity conservation rules constrain the dynamics. 
In particular, they prevent a single QP to move in the presence of such a perturbation. 
Consequently the 1-QP energy level of the toric code {9 — 0) does not give rise to 
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dispersive bands but is simply renormalized when h± 0. The corresponding gap 
A is shown in figure 10 (upper panel) and does not vanish at the self-dual point. 
However, this observation is compatible with the existence of a first-order transition 
that is due to level crossings which cannot be captured by analyzing low-energy levels. 
This situation is exactly the same as the one discussed in [18] for the Z2 toric code in 
a transverse field. As in [18], 2-QP bound states are either pinned or mobile in one 
dimension only, while the simplest two-dimensional dispersing object is a 4-QP bound 
state. 
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Figure 10. Lower panel : ground-state energy per site eo of H{0, h±,0) as a 
function of 8 (with J = ^ cos 9 and h± = | sin 6). Inset : zoom of the kink at 
6 = it/A indicating the presence of a first-order phase transition. Upper panel ; 
1-QP gap A as a function of 6 computed by pCUT. At order 7, other dlog-Pade 
approximants of A are very close to the one shown here. 



5. Conclusion and perspectives 

In this work, we have introduced exactly solvable models with Zjv (Abelian) anyons 
that generalize Kitaev's famous toric code [ ] or, equivalently, Wen's plaquette model 
[27]. Our main motivation was to probe the robustness of 1jn>2 topologically-ordered 
phases following recent works on the case N = 2 [14-:H ]. This question is crucial 
since such phases are believed to be protected against external perturbations up to 
an order proportional to the typical system size. Besides, this robustness has been 
proven recently for any local perturbation [4, ■",]. However, as early noticed by Kitaev 
in his seminal paper [ ] : "Of course, the perturbation should be small enough, or else 
a phase transition may occur." 

To investigate the limits of the standard aforementioned pcrturbative argument, 
we added local perturbations to this Z jv toric code and we obtained several important 
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results. First of all, for specific choices of the perturbation, the perturbed Zjv 
Hamiltonian can be mapped onto the two-dimensional iV-state quantum clock model in 
a transverse field. This mapping generalizes the correspondence between the Z2 toric 
code in a parallel field and the transverse-field Ising model [f 4, 1-^] to arbitrary N. For 
TV = 3 and antiferromagnetic couplings we found evidence for a second-order phase 
transition but we have not been able to determine accurately its universality class. Let 
us simply note that our estimate of the critical exponent zv is compatible with that 
of the three-dimensional classical XY model describing the three-dimensional three- 
state classical antiferromagnetic Potts model (see, e.g., [57]). Second, we have also 
shown that for iV = 3, multi-critical points are likely present in the phase diagram of 
H{hx,0, hz) [defined in (22)] although, here again, it would be important to analyze 
them with complementary tools. In particular, numerical simulations would be as 
valuable as for = 2 [ ]. In addition, as already observed for the N — 2 case 
[18], we have shown that self-dual Hamiltonians may arise provided the perturbating 
operators satisfy some constraints that we derived for arbitrary TV. In this case, the 
self-duality is associated with a dimensional reduction and, for iV = 3, we found that 
a first-order transition occurs at the self-dual point (as was also the case for = 2). 

From a methodological point of view, the combination of high-order perturbation 
theory (pCUT) and variational techniques (iPEPS) has been shown to be very efficient 
in a domain where few alternative approaches exist. From that respect, let us mention 
that improvement of the variational ansatz could certainly be achieved by taking into 
account the gauge symmetry in the tensor network as recently discussed in [20, 58, 59]. 

To conclude, we would like to mention related problems that would be worth 
investigating in order to deepen our understanding of topological phases' fragility. In 
the simplest case = 2, it would already be worthwhile considering a perturbed 
toric code on a lattice which is not self-dual such as, for instance, the honeycomb 
lattice. Indeed, in this case, the role of charges and fluxes degrees of freedom cannot 
be exchanged by a simple lattice transformation and it is likely that such a model could 
lead to a non-trivial phase diagram in the presence of a uniform magnetic field. One 
may also directly perturb Wen's plaquette model with arbitrary local perturbations. 
Indeed, as shown in [(>()] for N = 2, although Wen's and Kitaev's models are (almost) 
the same in the absence of perturbation, they display very different properties in the 
presence of a uniform magnetic field. To go beyond iV = 2, we believe that the study 
of the breakdown of Z^v topological phases initiated here would greatly benefit from a 
large- A'' analysis that might be performed in a field-theoretical framework. Although 
it is clearly beyond the scope of this paper it may shed light on several issues as 
for instance the universality class that can be met at the boundaries of topological 
phases. Let us mention that recent studies suggested the possibility that conformal 
quantum critical points may describe these phase transitions (see for instance [Gl, 02]). 
Finally, the most important step certainly consists in analyzing the robustness of non- 
Abelian topological phases [22, 63] which are of direct interest for topological quantum 
computation but, undoubtedly, these quantum objects are much more tricky to handle. 
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Appendix A. Series expansion in the limit \hz \ ^ J for hx = h± = 



Setting J = 1/3, the series expansion at order 11 of the ground-state energy per site 
of H{0,0,hz) reads 
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3 z z 2 ^ 36 ^ 144 ^ 27 ^ 311040 ^ 
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The series expansion at order 11 of the 1-QP gap reads 
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Table Al. Tabic of critical values and gap exponents obtained from dlog- 
Padc [n, m] extrapolation of the 1-QP gap A. Defective approxiniants, i.e., the 
ones with spurious poles between zero and /j^, arc marked with an asterisk. For 
/i^ > we do not give the corresponding exponent since, in this case, a first-order 
transition occurs at h* < /i^. 
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Appendix B. Series expansion in the limit |/ijf |, j/i^j <^ J for h± = 

Setting J = 1/3, the series expansion at order 7 of the ground-state energy per site of 
H{hx,0,hz) reads 

2 17 3 847 21 18407 

eo= -^-^S2-Ss--S, + -P,-—S, + -P,S,-^^Se 

33 933 15290 1477451 72039 

+ 32^^+ 40"^'^^^"^'^^+^9200"^'^'^^ + ^^^^^' ^^-^^ 

with Sn = h'^+ /i| and P„ = h'^^^hz^^. 

For iV = 3 and for nonvanishing hx and /i^, the 1-QP gap is defined as 
A = inin(Ac, Aj) where Ac (Aj^) denotes the charge (flux) gap. The series expansion 
at order 7 for the charge gap reads 
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The flux gap A/ is simply obtained from Ac by exchanging hx and /iz- 

Appendix C. Series expansion in the limit < <C J for hx = hz = 

Setting J = gcos0, h± = |sin0, and denoting t = ta,n0 = h±/J <C 1, the series 
expansion at order 7 of the ground-state energy per site of H{0, h±,0) reads 
eo 2 1 ^2 1 .3 599 ^4 209 ^- 896417 



cos6» 3 18 216 272160 259200 2204496000 

7184765443 ^7 (C.l) 



33331979520000 
The series expansion at order 7 of the 1-QP gap reads 
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Since this model is self-dual, the corresponding quantities in the opposite limit 
h± ^ J > are obtained by changing 6 into f — 0- 
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